##############################################################################
# File-Name: 03_analysis_sm_additional_results.r
# Purpose: Analysis primary outcomes
# Machine: macOS Ventura 13.5.1
# R version 4.3.1 (2023-06-16)
##############################################################################

# Packages ----------------------------------------------------------------
pacman::p_load(tidyverse, scales,  knitr, srvyr, extrafont, 
               rebus, broom, tidyr, lubridate, here, ggdist, ggtext, ggalt, janitor, 
               broom, tidyr,  stargazer, janitor, estimatr, list, marginaleffects, 
               modelsummary, kableExtra)


# Source Auxiliary Code ---------------------------------------------------
source(here("code","_utils.r"))

# Data wrangling ----------------------------------------------------------

# open
d <- readRDS(here("data","all_processed_data.rds"))

# create other datasets for partisans of PT and anti-partisans of PT
pt <- d %>% filter(partisans_pt==1)
antipt <- d %>% filter(anti_partisans_pt==1)

# for voting preference
# Leftists voters
dleft <- d %>% filter(w1_vote_first_recode=="Left")
dr <- d %>% filter(w1_vote_first_recode=="Right")
dm <- d %>% filter(str_detect(w1_vote_first_recode, "Others"))


# Fig 5 -----

# Recall of false items
list_false_items <- c("false_1_exp", "false_2_exp", "false_3_exp", "false_4_exp")
list_itt <- map(list_false_items, ~ func_ols_base(.x, data=d))
list_ittc <- map(list_false_items, ~ func_ols_covariates(.x, data=d))
list_cace <- map(list_false_items, ~ func_cace(.x, data=d))

# plot all together
output = flatten(list(list_itt, list_ittc, list_cace))
models= c(rep("ITT", 4), rep("Cov-ITT", 4), rep("CACE", 4))
dep_var_label= rep(c("False Item 1", "False Item 2","False Item 3", "False Item 4"), 3)

# plots
false_exp = plot_model(output,
                       models, 
                       dep_var_label, 
                       title="", 
                       caption= "Standardized Treatment Effects on exposure to False headlines by items") 
## Recall of true items
list_true_items <- c("true_1_exp", "true_2_exp", "true_3_exp", "true_4_exp")
list_itt <- map(list_true_items, ~ func_ols_base(.x, data=d))
list_ittc <- map(list_true_items, ~ func_ols_covariates(.x, data=d))
list_cace <- map(list_true_items, ~ func_cace(.x, data=d))

# plot all together
output = flatten(list(list_itt, list_ittc, list_cace))
models= c(rep("ITT", 4), rep("Cov-ITT", 4), rep("CACE", 4))
dep_var_label= rep(c("True Item 1", "True Item 2","True Item 3", "True Item 4"), 3)

# plots
true_exp <- plot_model(output,
                       models, 
                       dep_var_label, 
                       title="", 
                       caption= "Standardized Treatment Effects on exposure to False headlines by items") 

# get the data from the two models
res <- bind_rows(false_exp$data, true_exp$data)  %>% 
    mutate(dv=fct_rev(fct_inorder(dv)), 
           model=fct_rev(model))

# plot with the data combined
ggplot(res, 
       aes(y=dv, 
           x=estimate, 
           #color=party,
           #group=party,
           label=round(estimate, 2))) +
    geom_errorbar(aes(xmin=lb90, 
                      xmax=up90), 
                  size=2, width=0, alpha=.8, 
                  position=position_dodge(width = .8)) +
    geom_errorbar(aes(xmin=lb95, 
                      xmax=up95),  
                  size=1, width=.1, alpha=.5, 
                  position=position_dodge(width = .8)) +
    geom_point(size=12, shape=21, fill="white",
               position=position_dodge(width = .8), 
               stroke=2, alpha=1)  +
    geom_text(
        fontface = "italic", 
        color = "black",
        size=4, 
        position = position_dodge(width=.8)) +
    geom_vline(aes(xintercept=0), linetype="dashed", color="tomato2", alpha=.8) +
    ylab("") +
    xlab("") +
    labs(title="", 
         subtitle="", 
         caption="Standardized Treatment Effects on Belief Accuracy") +
    facet_grid(~model) +
    theme(legend.position = "bottom")



# save
save_function("sm_fig5.png")
save_function("sm_fig5.tiff")

# Fig 6 -----
## Separate effects of exposure by voting preferences -----

# DV: Exposure to false news
depvar = "false_news_exp"
false_true_itt <- map(list(dleft, dr, dm), ~ func_ols_base(depvar, data=.x))
false_true_ittc <- map(list(dleft, dr, dm),~ func_ols_covariates(depvar,.x))
false_true_cace <- map(list(dleft, dr, dm), ~ func_cace(depvar, data=.x))

# plot all together
output= flatten(list(false_true_itt, false_true_ittc, false_true_cace))
models= c(rep("ITT", 3), rep("Cov-ITT", 3), rep("CACE", 3))
dep_var_label=rep(c("Leftists", "Conservatives","Independents"), 3)

# plot
plot_model(output, models, dep_var_label,
           title="", 
           caption= "Standardized Treatment Effects on exposure to False headlines by partisan preferences") 

save_function("sm_fig6.png")
save_function("sm_fig6.tiff")


# Fig 7 ----
## Pro and Counter Attitudinal exposure  ---------------------------------------------

# Recoding exposure conditional on pro/counter attitudinal
d  <- d %>%
    mutate(pro_attitudinal_right=false_1_exp + false_3_exp, 
           counter_attitudinal_right=false_4_exp, 
           pro_attitudinal_left=false_4_exp, 
           counter_attitudinal_left=false_1_exp + false_3_exp)


# Leftists voters
dleft <- d %>% filter(w1_vote_first_recode=="Left")
list_false_items <- list("pro_attitudinal_left", "counter_attitudinal_left")
list_itt <- map(list_false_items, ~ func_ols_base(.x, data=dleft))
list_ittc <- map(list_false_items, ~ func_ols_covariates(.x, data=dleft))
list_cace <- map(list_false_items, ~ func_cace(.x, data=dleft))

# plot all together
output = flatten(list(list_itt, list_ittc, list_cace))
models= c(rep("ITT", 2), rep("Cov-ITT", 2), rep("CACE", 2))
dep_var_label= rep(c("Pro Attitudinal", "Counter Attitudinal"), 3)

# plots
left <- plot_model(output,
                   models, 
                   dep_var_label, 
                   title="A", 
                   caption="Standardized Treatment Effects on exposure to False headlines by partisan preferences") 

# Conservative voters
dr <- d %>% filter(w1_vote_first_recode=="Right")
list_false_items <- list("pro_attitudinal_right", "counter_attitudinal_right")
list_itt <- map(list_false_items, ~ func_ols_base(.x, data=dr))
list_ittc <- map(list_false_items, ~ func_ols_covariates(.x, data=dr))
list_cace <- map(list_false_items, ~ func_cace(.x, data=dr))

# plot all together
output = flatten(list(list_itt, list_ittc, list_cace))
models= c(rep("ITT", 2), rep("Cov-ITT", 2), rep("CACE", 2))
dep_var_label= rep(c("Pro Attitudinal", "Counter Attitudinal"), 3)

# plots
right <- plot_model(output,
                    models, 
                    dep_var_label, 
                    title="", 
                    caption="Standardized Treatment Effects on exposure to False headlines by partisan preferences") 

# combine both graphs
## colors
leftc=wesanderson::wes_palette("BottleRocket2")[[2]]
rightc=wesanderson::wes_palette("BottleRocket2")[[1]]

# plot
bind_rows(left$data %>% mutate(party="Leftists"), 
          right$data %>% mutate(party="Conservatives")) %>%
    ggplot(aes(y=model, 
               x=estimate, 
               color=party,
               group=party,
               label=round(estimate, 2))) +
    geom_errorbar(aes(xmin=lb90, 
                      xmax=up90), 
                  size=2, width=0, alpha=.8, 
                  position=position_dodge(width = .8)) +
    geom_errorbar(aes(xmin=lb95, 
                      xmax=up95),  
                  size=1, width=.1, alpha=.5, 
                  position=position_dodge(width = .8)) +
    geom_point(size=12, shape=21, fill="white",
               position=position_dodge(width = .8), 
               stroke=2, alpha=1)  +
    geom_text(
        fontface = "italic", 
        color = "black",
        size=4, 
        position = position_dodge(width=.8)) +
    geom_vline(aes(xintercept=0), linetype="dashed", color="tomato2", alpha=.8) +
    ylab("") +
    xlab("") +
    labs(title="", 
         subtitle="", 
         caption="Standardized Treatment Effects on exposure to False headlines by partisan preferences") +
    scale_color_manual(values=c(rightc, leftc),  name="Partisan Preferences", 
                       guide = guide_legend(reverse = TRUE, 
                                            override.aes = list(size=3), 
                                            title.position="top")) +
    facet_grid(~dv) +
    theme(legend.position = "bottom")


save_function("sm_fig7.png")
save_function("sm_fig7.tiff")



# Fig 8  ------------------------------------------------------------------
## Truth discernment as DV------------
# Truth Discernment is the sum of accuracy for false and true
d$td <- d$false_false_sum+d$true_true_sum
depvar = "td"

# models
td_itt <- func_ols_base(depvar, data=d)
td_ittc <- func_ols_covariates(depvar, data=d)
td_cace <- func_cace(depvar, data=d)

# plot all together
output = list(td_itt, td_ittc, td_cace)
models = list("ITT", "Cov-ITT", "CACE")
dep_var_label = c(rep("Truth Discernment", 3))

# plot
plot_model(output, models, dep_var_label, 
           title="", 
           caption= "Standardized Treatment Effects on Truth Discerniment")


save_function("sm_fig8.png")
save_function("sm_fig8.tiff")

# Fig9 ----
## Models with results by item  -------

# accuracy on false items
list_false_items <- c("false_1_false", "false_2_false", "false_3_false", "false_4_false")
list_itt <- map(list_false_items, ~ func_ols_base(.x, data=d))
list_ittc <- map(list_false_items, ~ func_ols_covariates(.x, data=d))
list_cace <- map(list_false_items, ~ func_cace(.x, data=d))

# plot all together
output = flatten(list(list_itt, list_ittc, list_cace))
models= c(rep("ITT", 4), rep("Cov-ITT", 4), rep("CACE", 4))
dep_var_label= rep(c("False Item 1", "False Item 2","False Item 3", "False Item 4"), 3)

# plots
false_accuracy <- plot_model(output,
                             models, 
                             dep_var_label, 
                             title="", 
                             caption= "Standardized Treatment Effects on Beliefs Accuraccy by False Item")


# accuracy on true items
list_true_items <- c("true_1_true", "true_2_true", "true_3_true", "true_4_true")
list_itt <- map(list_true_items, ~ func_ols_base(.x, data=d))
list_ittc <- map(list_true_items, ~ func_ols_covariates(.x, data=d))
list_cace <- map(list_true_items, ~ func_cace(.x, data=d))

# plot all together
output = flatten(list(list_itt, list_ittc, list_cace))
models= c(rep("ITT", 4), rep("Cov-ITT", 4), rep("CACE", 4))
dep_var_label= rep(c("True Item 1", "True Item 2","True Item 3", "True Item 4"), 3)

# plots
true_accuracy <- plot_model(output,
                            models, 
                            dep_var_label, 
                            title="", 
                            caption= "Standardized Treatment Effects on Beliefs Accuraccy by True Items") 

# combine datasets
res <- bind_rows(true_accuracy$data, false_accuracy$data)  %>% 
    mutate(dv=fct_rev(fct_inorder(dv)), 
           model=fct_rev(model))

# plot
ggplot(res, 
       aes(y=dv, 
           x=estimate, 
           #color=party,
           #group=party,
           label=round(estimate, 2))) +
    geom_errorbar(aes(xmin=lb90, 
                      xmax=up90), 
                  size=2, width=0, alpha=.8, 
                  position=position_dodge(width = .8)) +
    geom_errorbar(aes(xmin=lb95, 
                      xmax=up95),  
                  size=1, width=.1, alpha=.5, 
                  position=position_dodge(width = .8)) +
    geom_point(size=12, shape=21, fill="white",
               position=position_dodge(width = .8), 
               stroke=2, alpha=1)  +
    geom_text(
        fontface = "italic", 
        color = "black",
        size=4, 
        position = position_dodge(width=.8)) +
    geom_vline(aes(xintercept=0), linetype="dashed", color="tomato2", alpha=.8) +
    ylab("") +
    xlab("") +
    labs(title="", 
         subtitle="", 
         caption="Standardized Treatment Effects on Belief Accuracy") +
    facet_grid(~model) +
    theme(legend.position = "bottom")

save_function("sm_fig9.png")
save_function("sm_fig9.tiff")


# Fig 10 ------------------------------------------------------------------

# Dependent Variable: Sum of false news accuracy statements
depvar = "false_false_sum"
false_true_itt <- map(list(dleft, dr, dm), ~ func_ols_base(depvar, data=.x))
false_true_ittc <- map(list(dleft, dr, dm),~ func_ols_covariates(depvar,.x))
false_true_cace <- map(list(dleft, dr, dm), ~ func_cace(depvar, data=.x))

# plot all together
output= flatten(list(false_true_itt, false_true_ittc, false_true_cace))
models= c(rep("ITT", 3), rep("Cov-ITT", 3), rep("CACE", 3))
dep_var_label=rep(c("Leftists", "Conservatives","Independents"), 3)

# plot
plot_model(output, models, dep_var_label, 
           title="A: False Rumors Accuracy", 
           caption= "Standardized Treatment Effects on False rumors accuracy conditional on partisan preferences")

save_function("sm_fig10a.png")
save_function("sm_fig10a.tiff")

# Dependent Variable: Sum of true news accuracy statements
depvar = "true_true_sum"
func_t.test(depvar, data=dleft)
false_true_itt <- map(list(dleft, dr, dm), ~ func_ols_base(depvar, data=.x))
false_true_ittc <- map(list(dleft, dr, dm),~ func_ols_covariates(depvar,.x))
false_true_cace <- map(list(dleft, dr, dm), ~ func_cace(depvar, data=.x))

# plot all together
output= flatten(list(false_true_itt, false_true_ittc, false_true_cace))
models= c(rep("ITT", 3), rep("Cov-ITT", 3), rep("CACE", 3))
dep_var_label=rep(c("Leftists", "Conservatives","Independents"), 3)

# plot
plot_model(output, models, dep_var_label, 
           title="B: True News Accuracy", 
           caption= "Standardized Treatment Effects on True news accuracy conditional on partisan preferences")

save_function("sm_fig10b.png")
save_function("sm_fig10b.tiff")


# Fig 11 and Fig 12 ------------------------------------------------------------------
## Polarization Outcomes
# All
list_politems_items <- list("pol_all_zcore", 
                            "distance_candidates",
                            "affective_pol_diff_all", 
                            "social_pol_agg", 
                            "mean_abs_policy_polarization")

list_pol_itt <- map(list_politems_items, ~ func_ols_base(.x, data=d))
list_pol_ittc <- map(list_politems_items, ~ func_ols_covariates(.x, data=d))
list_pol_cace <- map(list_politems_items, ~ func_cace(.x, data=d))

# plot all together
output = flatten(list(list_pol_itt, list_pol_ittc, list_pol_cace))
models= c(rep("ITT", 5), rep("Cov-ITT", 5), rep("CACE", 5))
dep_var_label= rep(c("Polarization \n Index", 
                     "Ideological \n Polarization",
                     "Affective \n Polarization", 
                     "Social \n Polarization", 
                     "Issue \n Polarization"), 3)


# subjective well-being - All items

# double check the index
list_wb_items <- list("swb_zcore", # zscore index
                      "q_well_being_happy_numeric"  ,
                      "q_well_being_depressed_numeric",
                      "q_well_being_anxious_numeric",
                      "q_well_being_isolated_numeric", 
                      "q_well_being_satisfied_numeric" )

# run the models
list_swb_itt <- map(list_wb_items, ~ func_ols_base(.x, data=d))
list_swb_ittc <- map(list_wb_items, ~ func_ols_covariates(.x, data=d))
list_swb_cace <- map(list_wb_items, ~ func_cace(.x, data=d))



### Combine Polarization outcomes and Subjective well-being -------
list_bind_itt <- c(list_pol_itt, list(list_swb_itt[[1]]))
list_bind_ittc <- c(list_pol_ittc, list(list_swb_ittc[[1]]))
list_bind_cace <- c(list_pol_cace, list(list_swb_cace[[1]]))

# plot all together
output = flatten(list(list_bind_itt, list_bind_ittc,list_bind_cace))
models= c(rep("ITT", 6), rep("Cov-ITT", 6), rep("CACE", 6))
dep_var_label= rep(c("Polarization Index", 
                     "Ideological \n Polarization",
                     "Affective Polarization", 
                     "Social Polarization", 
                     "Issue Polarization", 
                     "Subjective Well-Being Index"), 3)

# only index
index <- plot_model(output,
                    models, 
                    dep_var_label, 
                    title="", 
                    caption ="") 

# all polarization
pol<- plot_model(output,
                 models, 
                 dep_var_label, 
                 title="", 
                 caption ="") 


# plot with each of the outcomes
pol$data <- pol$data %>%
    filter(!str_detect(dv, "Index"))

pol +
    labs(caption="Standardized Treatment Effects on Polarization Items", 
         title="")

save_function("sm_fig11.png")
save_function("sm_fig11.tiff")

# swb outcomes -- Not using them recoded because the original scale is more intuitive
d <- d %>% 
    mutate_at(vars(c(q_well_being_anxious, 
                     q_well_being_depressed, 
                     q_well_being_isolated, 
                     q_well_being_happy, 
                     q_well_being_satisfied)), as.numeric)

list_wb_items <- list("q_well_being_happy"  ,
                      "q_well_being_satisfied" ,
                      "q_well_being_depressed",
                      "q_well_being_anxious",
                      "q_well_being_isolated")

# get only partisans
list_swb_itt <- map(list_wb_items, ~ func_ols_base(.x, data=d))
list_swb_ittc <- map(list_wb_items, ~ func_ols_covariates(.x, data=d))
list_swb_cace <- map(list_wb_items, ~ func_cace(.x, data=d))

# plot all together
output = flatten(list(list_swb_itt, list_swb_ittc, list_swb_cace))
models= c(rep("ITT", 5), rep("Cov-ITT", 5), rep("CACE", 5))
dep_var_label= rep(c("Happy",
                     "Satisfied",
                     "Depressed", 
                     "Anxious", 
                     "Isolated"), 3)

# plots
swb <- plot_model(output,
                  models, 
                  dep_var_label, 
                  title="", 
                  subtitle= "", 
                  caption="Treatment Effects on Subjective Well-Being items") 

# rename the variables
swb$data <- swb$data %>%
    mutate(dv=fct_relevel(dv, "Happy", "Satisfied", "Depressed", "Anxious", "Isolated")) 

swb

save_function("sm_fig12.png")
save_function("sm_fig12.tiff")

# Figure 15 ---------------------------------------------------------------

# Combining the counts of exposure per participant
d_2 <- d %>% select(response_id, exp, false_news_exp, true_news_e) %>%
    pivot_longer(cols=-c(response_id,exp), 
                 names_to="headline", 
                 values_to="dv") %>%
    mutate(mod=ifelse(str_detect(headline, "false"), "false", "true"), 
           mod=fct_rev(mod))

# model the interactive effect
model2 = lm_robust(dv ~ exp*mod, data=d_2, 
                   clusters = response_id)  


# estimate marginal effect
res2 <- avg_slopes(
    model2,
    variables = "exp",
    by = c("mod")) %>% tidy() %>%
    mutate(up90=estimate + 1.64*std.error,
           up95=estimate+ 1.96*std.error, 
           lb90=estimate - 1.64*std.error, 
           lb95=estimate - 1.96*std.error, 
           moderator_f=fct_rev(ifelse(mod=="true", "True News Familiarity",
                                      "False Rumors Familiarity")))

# check pvalues
hypothesis_test <- car::linearHypothesis(model2, c("exptreatment+exptreatment:modfalse=exptreatment"))
pvalues <- hypothesis_test$`Pr(>Chisq)`[2]

# Plot Marginal Effects Interactive
ggplot(res2, 
       aes(y=estimate, 
           x=moderator_f)) +
    geom_errorbar(aes(ymin=lb90, 
                      ymax=up90), 
                  size=2, width=0, alpha=.8, 
                  position=position_dodge(width = .8), 
                  color=colorwes) +
    geom_errorbar(aes(ymin=lb95, 
                      ymax=up95),  
                  size=1, width=.1, alpha=.5, 
                  position=position_dodge(width = .8), 
                  color=colorwes) +
    geom_point(size=14, shape=21, fill="white",
               position=position_dodge(width = .6), 
               stroke=2, 
               color=colorwes, alpha=1) +
    labs(title="", 
         y="Marginal Effects", 
         x="", 
         caption="") +
    geom_segment(aes(x=1, xend=2, y=0, yend=0), 
                 size=.8, 
                 color="black") +
    geom_segment(aes(x=1, xend=1, y=-0, yend=-0.01), 
                 size=.8, 
                 color="black") +
    geom_segment(aes(x=2, xend=2, y=-0, yend=-0.01), 
                 size=.8, 
                 color="black") +
    geom_text(aes(x=1.5, y=0.02, 
                  label=paste0("p-value=", round(pvalues, digits=3)), family=""),
              size=4, family=my_font, color = "black")  +
    geom_text(aes(label=round(estimate, digits=2)), 
              fontface = "italic", 
              fill="white", 
              color = "black",
              size=4, 
              position = position_dodge(.6))

save_function("sm_fig15.png")
save_function("sm_fig15.tiff")


# Figure 19 ---------------------------------------------------------------
# Comparison of different measure of polarization

# Polarization based on voting choice
dvote <- d %>% filter(q_runoff=="Left"|q_runoff=="Right")

# Polarization based on Feelings about the PT 
dpartisans <- d %>% filter(partisans_pt==1|anti_partisans_pt==1)

# Dependent variable: Polarization ------------------------------------------------------------

## All participants

# All
list_politems_items <- list("pol_all_zcore", 
                            "distance_candidates",
                            "affective_pol_diff_all", 
                            "social_pol_agg", 
                            "mean_abs_policy_polarization")

list_pol_itt <- map(list_politems_items, ~ func_ols_base(.x, data=d))
list_pol_ittc <- map(list_politems_items, ~ func_ols_covariates(.x, data=d))
list_pol_cace <- map(list_politems_items, ~ func_cace(.x, data=d))

# plot all together
output = flatten(list(list_pol_itt, list_pol_ittc, list_pol_cace))
models= c(rep("ITT", 5), rep("Cov-ITT", 5), rep("CACE", 5))
dep_var_label= rep(c("Polarization \n Index", 
                     "Ideological \n Polarization",
                     "Affective \n Polarization", 
                     "Social \n Polarization", 
                     "Issue \n Polarization"), 3)


# plots
s <- plot_model(output,
                models, 
                dep_var_label, 
                title="", 
                subtitle= "", 
                caption="Standardized Treatment Effects on Polarization") +
    theme(axis.text.x  = element_text(size=10), 
          strip.text = element_text(family = my_font, 
                                    color = "#22211d",
                                    size = 10, face="bold")
    )


## Models with only partisans ----------------------------------------------
list_pol_itt <- map(list_politems_items, ~ func_ols_base(.x, data=dvote))
list_pol_ittc <- map(list_politems_items, ~ func_ols_covariates(.x, data=dvote))
list_pol_cace <- map(list_politems_items, ~ func_cace(.x, data=dvote))

# plot all together
output = flatten(list(list_pol_itt, list_pol_ittc, list_pol_cace))
models= c(rep("ITT", 5), rep("Cov-ITT", 5), rep("CACE", 5))
dep_var_label= rep(c("Polarization \n Index", 
                     "Ideological \n Polarization",
                     "Affective \n Polarization", 
                     "Social \n Polarization", 
                     "Issue \n Polarization"), 3)


# plots
s_vote <- plot_model(output,
                     models, 
                     dep_var_label, 
                     title="", 
                     subtitle= "", 
                     caption="Standardized Treatment Effects on Polarization") +
    theme(axis.text.x  = element_text(size=10), 
          strip.text = element_text(family = my_font, 
                                    color = "#22211d",
                                    size = 10, face="bold")
    )


## Only with Feeling about Petistas ----------------------------------------------

list_pol_itt <- map(list_politems_items, ~ func_ols_base(.x, data=dpartisans))
list_pol_ittc <- map(list_politems_items, ~ func_ols_covariates(.x, data=dpartisans))
list_pol_cace <- map(list_politems_items, ~ func_cace(.x, data=dpartisans))


# plot all together
output = flatten(list(list_pol_itt, list_pol_ittc, list_pol_cace))
models = c(rep("ITT", 5), rep("Cov-ITT", 5), rep("CACE", 5))
dep_var_label= rep(c("Polarization \n Index", 
                     "Ideological \n Polarization",
                     "Affective \n Polarization", 
                     "Social \n Polarization", 
                     "Issue \n Polarization"), 3)


# plots
s_partisans <- plot_model(output,
                          models, 
                          dep_var_label, 
                          title="", 
                          subtitle= "", 
                          caption="Standardized Treatment Effects on Polarization") +
    theme(axis.text.x  = element_text(size=10), 
          strip.text = element_text(family = my_font, 
                                    color = "#22211d",
                                    size = 10, face="bold")
    )



# Combine plots -----------------------------------------------------------

# Edit
s$data <- s$data %>%
    mutate(dv=fct_relevel(dv, "Polarization \n Index", 
                          "Ideological \n Polarization",
                          "Affective \n Polarization", 
                          "Social \n Polarization", 
                          "Issue \n Polarization"),
           p="Full Sample") 

s_vote$data <- s_vote$data %>%
    mutate(dv=fct_relevel(dv, "Polarization \n Index", 
                          "Ideological \n Polarization",
                          "Affective \n Polarization", 
                          "Social \n Polarization", 
                          "Issue \n Polarization"),
           p="Partisan Voters") 

s_partisans$data <- s_partisans$data %>%
    mutate(dv=fct_relevel(dv, "Polarization \n Index", 
                          "Ideological \n Polarization",
                          "Affective \n Polarization", 
                          "Social \n Polarization", 
                          "Issue \n Polarization"),
           p="Positive and Negative \n Workers Party") 



bind_rows(s$data, s_vote$data, s_partisans$data) %>%
    ggplot(aes(y=model, 
               x=estimate, 
               color=p,
               group=p,
               label=round(estimate, 2))) +
    geom_errorbar(aes(xmin=lb90, 
                      xmax=up90), 
                  size=2, width=0, alpha=.8, 
                  position=position_dodge(width = .8)) +
    geom_errorbar(aes(xmin=lb95, 
                      xmax=up95),  
                  size=1, width=.1, alpha=1, 
                  position=position_dodge(width = .8)) +
    geom_point(size=2, shape=21, fill="white",
               position=position_dodge(width = .8), 
               stroke=2, alpha=1)  +
    # geom_text(
    #   fontface = "italic", 
    #   color = "black",
    #   size=4, 
    #   position = position_dodge(width=.8)) +
    geom_vline(aes(xintercept=0), linetype="dashed", color="tomato2", alpha=.8) +
    ylab("") +
    xlab("Standardized Treatment Effects") +
    labs(title="", 
         subtitle="", 
         caption="") +
    scale_color_manual(values=rev(c("#969696", "#525252", "#000000")), 
                       name="Samples:", 
                       guide = guide_legend(reverse = TRUE, 
                                            override.aes = list(size=3), 
                                            title.position="left")) +
    facet_grid(~dv) +
    theme(legend.position = "bottom", 
          axis.title.x = element_text(size=12))

save_function("sm_fig19.png")
save_function("sm_fig19.tiff")



# Regression Tables --------------------------------------------------------

## table 9 - Dependent variable: Self-reported exposure --------------------------------------------------

# rotate the variable so that higher values means more exposure
d <- d %>%
    mutate(q_fake_news=fct_rev(q_fakenews), # relevel
           q_fake_news_no=as.numeric(q_fake_news)) # convert to numeric

#Select dependent variable
depvar <- "q_fake_news_no"

# Estimate the models
self_e_itt <- func_ols_base(depvar, data=d, model=TRUE)
self_e_ittc <- func_ols_covariates(depvar,d, model=TRUE)
self_e_cace <- func_cace(depvar, data=d, model=TRUE)

# Model with correct interaction terms
new_names = c ("(Intercept)"  = "Intercept",  
               "exptreatment" ="Treatment",
               "compliance_one_side"  = "Treatment")


# options for modelsummary
options("modelsummary_format_numeric_latex" = "plain")


# Save
library(modelsummary)
## latex version
# modelsummary(list("ITT"=self_e_itt, 
#                   "Cov-ITT"=self_e_ittc, 
#                   "CACE" = self_e_cace),
#              # if a latex issue is raised, just comment out everything below, and run only the modelsummary to see the results
#              output="latex",
#              fmt = 3,
#              stars = TRUE,  
#              coef_map = new_names, 
#              title = "Regression Models: Deactivation Treatment Effects on Self-Reported Exposure to Misinformation on WhatsApp. The outcome variable is a nominal scale measuring self-reported exposure to misinformation. Results indicate robust negative descrease on self-reported exposure. ",
#              gof_omit = 'AIC|BIC') %>%
#     footnote(general = "Robust standard errors in Parentheses", 
#              threeparttable = TRUE)  %>%
#     save_kable(file = "output/sm_tab9.tex")

## if a latex issue is raised, you can check results with this: 
modelsummary(list("ITT"=self_e_itt, 
                  "Cov-ITT"=self_e_ittc, 
                  "CACE" = self_e_cace), 
             fmt = 3,
             stars = TRUE,  
             coef_map = new_names, 
             title = "Regression Models: Deactivation Treatment Effects on Self-Reported Exposure to Misinformation on WhatsApp. The outcome variable is a nominal scale measuring self-reported exposure to misinformation. Results indicate robust negative descrease on self-reported exposure. ",
             gof_omit = 'AIC|BIC') 

## Trust -----

# Table 13 - Regression table for Trust outcomes -------------------------------------------------------------------

# get all trust items
list_trust_items <- list("trust_zcore", 
                         "q_trust_government"  ,
                         "q_trust_congress" ,
                         "q_trust_electoral_authorities",
                         "q_trust_globo",
                         "q_trust_news_channels")

# models
list_trust_itt <- map(list_trust_items, ~ func_ols_base(.x, data=d))
list_trust_ittc <- map(list_trust_items, ~ func_ols_covariates(.x, data=d, model=TRUE))
list_trust_cace <- map(list_trust_items, ~ func_cace(.x, data=d))

# Model with correct interaction terms
new_names = c ("(Intercept)"  = "Intercept",  
               "exptreatment" = "Treatment")

# table
# modelsummary(list("Trust Score"=list_trust_ittc[[1]], 
#                   "Trust Government"=list_trust_ittc[[2]], 
#                   "Trust Congress" = list_trust_ittc[[3]], 
#                   "Trust Electoral Authorities"=list_trust_ittc[[4]],
#                   "Trust Globo (TV News Media)"=list_trust_ittc[[5]], 
#                   "Trust News Outlets"=list_trust_ittc[[6]]),
#              # if a latex issue is raised, just comment out everything below, and run only the modelsummary to see the results
#              output="latex",
#              fmt = 3,
#              stars = TRUE,  
#              coef_map = new_names, 
#              title = "Regression Models: Covariate Adjusted Models for Deactivation Treatment Effects on Trust Measures. The survey data uses a seven point nominal scale measuring trust in a range of institutions and news media.  The outcome for the models uses a standardized z-score for the indicator trust measures, and the sum of the z-scores for the index. Our results do not recover any statistically significant effect of the deactivation on standard measures for trust",
#              gof_omit = 'AIC|BIC') %>%
#     footnote(general = "Robust standard errors in Parentheses", 
#              threeparttable = TRUE)   %>%
#     save_kable(file = "output/sm_tab13.tex")

## if a latex issue is raised, you can check results with this: 
modelsummary(list("Trust Score"=list_trust_ittc[[1]], 
                  "Trust Government"=list_trust_ittc[[2]], 
                  "Trust Congress" = list_trust_ittc[[3]], 
                  "Trust Electoral Authorities"=list_trust_ittc[[4]],
                  "Trust Globo (TV News Media)"=list_trust_ittc[[5]], 
                  "Trust News Outlets"=list_trust_ittc[[6]]),
             # if a latex issue is raised, just comment out everything below, and run only the modelsummary to see the results
             fmt = 3,
             stars = TRUE,  
             coef_map = new_names, 
             title = "Regression Models: Covariate Adjusted Models for Deactivation Treatment Effects on Trust Measures. The survey data uses a seven point nominal scale measuring trust in a range of institutions and news media.  The outcome for the models uses a standardized z-score for the indicator trust measures, and the sum of the z-scores for the index. Our results do not recover any statistically significant effect of the deactivation on standard measures for trust",
             gof_omit = 'AIC|BIC')


# Table 14 - Regression table for Substitution effects ------------------------------------------------------------

# convert substitution outcomes to numeric
d <- d %>% 
    mutate_at(vars(contains("q_subs")), as.numeric) 

# get dvs names
list_trust_items <- list("q_subs_other_social_media", 
                         "q_subs_reading_news_internet"  ,
                         "q_subs_watching_tv" ,
                         "q_subs_offline_friends",
                         "q_subs_away_phone")

# Models
list_trust_itt <- map(list_trust_items, ~ func_ols_base(.x, data=d, model = TRUE))
list_trust_ittc <- map(list_trust_items, ~ func_ols_covariates(.x, data=d, model = TRUE))
list_trust_cace <- map(list_trust_items, ~ func_cace(.x, data=d))

# Model with correct interaction terms
new_names = c ("(Intercept)"  = "Intercept",  
               "exptreatment" = "Treatment")
# table - latex version
# modelsummary(list("Other Social Media"=list_trust_ittc[[1]], 
#                   "Reading News Online"=list_trust_ittc[[2]], 
#                   "Watching TV" = list_trust_ittc[[3]], 
#                   "Offline with Friends"=list_trust_ittc[[4]],
#                   "Away from my Phone"=list_trust_ittc[[5]]),
#              # if a latex issue is raised, just comment out everything below, and run only the modelsummary to see the results
#              output="latex",
#              fmt = 3,
#              stars = TRUE,  
#              coef_map = new_names, 
#              title = "Regression Models: Covariate Adjusted Models for Substitution effects of the Deactivation Treatment. The survey data uses a five point nominal scale asking how much less/much more time participants spent in a set of online and offline activities during the weeks of the experiment. Our results do not uncover an substitution of WhatsApp by other online activites. We find participants rely more on TV News as an consequence of the deactivation.",
#              gof_omit = 'AIC|BIC') %>%
#     footnote(general = "Robust standard errors in Parentheses", 
#              threeparttable = TRUE)   %>%
#     save_kable(file = "output/sm_tab14.tex")

## if a latex issue is raised, you can check results with this: 
modelsummary(list("Other Social Media"=list_trust_ittc[[1]], 
                  "Reading News Online"=list_trust_ittc[[2]], 
                  "Watching TV" = list_trust_ittc[[3]], 
                  "Offline with Friends"=list_trust_ittc[[4]],
                  "Away from my Phone"=list_trust_ittc[[5]]),
             fmt = 3,
             stars = TRUE,  
             coef_map = new_names, 
             title = "Regression Models: Covariate Adjusted Models for Substitution effects of the Deactivation Treatment. The survey data uses a five point nominal scale asking how much less/much more time participants spent in a set of online and offline activities during the weeks of the experiment. Our results do not uncover an substitution of WhatsApp by other online activites. We find participants rely more on TV News as an consequence of the deactivation.",
             gof_omit = 'AIC|BIC')
